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Abstract 

1 Motivation: 

The simulation of morphogenetic problems requires 
the simultaneous and coupled simulation of signalling 
and tissue dynamics. A cellular resolution of the 
tissue domain is important to adequately describe 
the impact of cell-based events, such as cell divi¬ 
sion, cell-cell interactions, and spatially restricted sig¬ 
nalling events. A tightly coupled cell-based mechano- 
regulatory simulation tool is therefore required. 

2 Results: 

We developed an open-source software framework for 
morphogenetic problems. The environment offers 
core functionalities for the tissue and signalling mod¬ 
els. In addition, the software offers great flexibility 
to add custom extensions and biologically motivated 
processes. Cells are represented as highly resolved, 
massless elastic polygons; the viscous properties of 
the tissue are modelled by a Newtonian fluid. The 
Immersed Boundary method is used to model the 
interaction between the viscous and elastic proper¬ 
ties of the cells, thus extending on the IBCell model. 
The fluid and signalling processes are solved using the 
Lattice Boltzmann method. As application examples 
we simulate signalling-dependent tissue dynamics. 

* simon. tanaka@bsse.ethz.ch 
t dagmar.iber@bsse.ethz.ch 


3 Availability: 

The documentation and source code are available 

on http://tanakas.bitbucket.org/lbibcell/ 
index.html 

4 Contact: 

dagmar. iber@bsse.ethz.ch 

5 Introduction 

During morphogenesis, tissue grows and self- 
organises into complex functional units such as or¬ 
gans. The process is tightly controlled, both by sig¬ 
nalling and by mechanical interactions. Long-range 
signalling interactions in the tissues can be medi¬ 
ated by diffusible substances, called morphogens, and 
by long-range cell processes ([35]). The dynamics of 
the diffusible factors can typically be well described 
by systems of continuous reaction-advection-diffusion 
partial differential equations (PDE). The appropriate 
tissue representation depends on the relevant time 
scale. For a homogeneous isotropic embryonic tis¬ 
sue, experiments show that the tissue is well approx¬ 
imated by a viscous fluid on long time scales (equi¬ 
libration after 30 minutes to several hours) and by 
an elastic material on short time scales (seconds to 
minutes) ([10]). However, biological control typically 
happens on a shorter time scale, and many cellular 
processes such as cell migration and adhesion, cell 
polarity, directed division, monolayer structures and 
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differentiation cannot be cast into a continuous for¬ 
mulation in a straight-forward way. A number of cell 
based simulation techniques at different scales and 
different level of detail have been developed to study 
these processes; here, we discuss main representatives 
for each category. 

The Cellular Potts model , introduced by m , is 
solved on a lattice, with each lattice point holding a 
generalized spin value denoting cell identity. Similar 
to the Ising model, Hamiltonian energy functions are 
formulated and minimized using a Metropolis algo¬ 
rithm. It has been applied to a multitude of prob¬ 
lems, and is implemented in the software Compu- 
CellSD ([39]). However, the correspondence between 
the biological problem and the Hamiltonian, the tem¬ 
perature and the time step is not always straightfor¬ 
ward. 

The subcellular element model divides cells into 
sub cellular elements, which are represented by com¬ 
putational particles. The elements interact via inter¬ 
acting potentials which are subject to modelling. The 
motion of the elements is governed by overdamped 
Langevin dynamics, such that the method is mesh- 
free. The framework was first introduced by [25] and 
later applied by [36] 137]. This approach allows for 
detailed biophysical modelling, both in 2D and 3D. 

The spheroid model developed by [7] assumes that 
cells in unstructured cell populations are similar to 
colloidal particles. The cells are modelled as point 
particles, hosting interaction potentials. Their mo¬ 
tion consist of a random and a directed movement. 
Neighboring cells form adhesive bonds, which are rep¬ 
resented using models borrowed from contact me¬ 
chanics, such as e.g. the Johnson-Kendall-Roberts 
model (ED- Many cellular processes such as cell 
shape change, division, death, lysis, cell-cell interac¬ 
tion and migration have been successfully translated 
into the spheroid model ([7]). Intra- and extracellu¬ 
lar diffusion has not yet been introduced and imple¬ 
mented. The spheroid model extends efficiently to 
3D, and it has been implemented in the open-source 
framework CellSys (USD- 

The vertex model uses polygons (or polyhedra in 
3D) to represent cells in densely packed tissues, e.g. 
in Drosophila wing disc epithelia (0). For each ver¬ 
tex, forces are computed - either via a potential or 


directly. The vertices are moved subsequently ac¬ 
cording to overdamped equations of motion, or via a 
Monte Carlo algorithm. The model is implemented 
in the open-source software Chaste am). 

The viscoelastic cell model (also called IBCell mod¬ 
els) presented in [34] and [33 uses the immersed 
boundary method ([29]) to represent individually de¬ 
formable cells as immersed elastic bodies. The cy¬ 
toplasm and the extracellular matrix and fluid are 
represented by a viscous incompressible fluid. In this 
framework, a vast amount of biological processes such 
as cell growth, cell division, apoptosis and polariza¬ 
tion has been realized. The model was applied to 
study tumor and epithelial dynamics. Due to the 
very high level of detail, the viscoelastic cell model 
is computationally expensive and has not yet been 
implemented in 3D. 

The software framework VirtualLeaf with explicit 
cell resolution, available in 2D, has been introduced 
in [22]. Although the cell representation is similar to 
vertex cell models, the dynamics is realized by min¬ 
imizing an Hamiltonian using a Monte Carlo algo¬ 
rithm. The model assumes rigid cell walls, which is 
appropriate for plant morphogenesis. 

For many morphogenetic phenomena, which arise 
from a tight interaction between the biomolecular sig¬ 
nalling and the tissue physics, an explicit computa¬ 
tional representation of the cell shapes is required. 
Here, we present a flexible software framework based 
on the IBCell model, which, as a novelty, permits 
to tightly couple biomolecular signalling models to a 
cell-resolved, physical tissue model. The core com¬ 
ponents and the general approach of the model are 
described in the second section. In the third sec¬ 
tion, the software and the main functionalities are 
described in detail. Application examples are given 
in the fourth chapter to demonstrate the framework’s 
capabilities. 

6 Approach 

Our approach permits the coupled simulation of tis¬ 
sue and signalling dynamics. To describe the tissue 
dynamics, the visco-elastic cell model needs to rep¬ 
resent both the cellular structures and their elastic 
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Figure 1: Algorithm Overview. The algorithm 
consists of three coupled layers. The geometry 
X (l,t) (top part, discussed in more detail in Figure 
[5]) is used to compute the forces F (/, t) acting on each 
of the geometry nodes. These forces, which do not 
necessarily coincide with a lattice point, are scattered 
to the fluid lattice (middle part) using the immersed 
boundary method kernel function, F (l,t) —> f (t). 
After advancing the fluid solver by one time step, the 
velocity is interpolated to the geometry node position 
using the same kernel function, u(x,t ) —)> L7 (/,£). 
The geometry nodes are moved according to their ve¬ 
locity U (l,t) and the iteration is restarted. The ve¬ 
locity u(x,t) of the fluid lattice is also copied to the 
reaction-advection-diffusion solvers (PDE), together 
with the position X (Z, t) —)> x (t) of the geometry. 
The state of the reaction-advection-diffusion solvers, 
which are used to model signalling, may be used to 
compute mass sources S (x,t) for the fluid solver. 


properties, as well as the viscous behaviour of the cy¬ 
toplasm and of the extracellular space surrounding 
the cells. The model therefore rests on three core 
parts: the representation of cells, the representation 
of the fluid and the fluid-structure interaction, and 
the coupling of the tissue part to the signalling model. 


To describe the interaction between the viscous fluid 
and the elastic structures, which are immersed in the 
fluid, we use the immersed boundary (IB) method 
([29]) as previously implemented in the visco-elastic 
cell model, also called IBcell model ([551 33]). To 
solve the viscous fluid behaviour, we use the Lattice 
Boltzmann method, which is an efficient mesoscopic 
numerical scheme, originally developed to solve fluid 
dynamics problems ([5]). The method has previously 
been successfully applied to reaction-diffusion equa¬ 
tions such as Turing systems ([32]), as well as to cou¬ 
pled scalar fields such as temperature ([13]). The 
method was for the first time combined with the im¬ 
mersed boundary method ([29]) by [9], and has later 
been used to study red blood cells in flow by [ 421 . 
In the following, we provide an overview of the im¬ 
plemented methods; the implementation details are 
given in section [7] 

6.1 Cell Representation 

Cells are represented as massless, purely elastic struc¬ 
tures, which are described by sets of geometry points 
forming polygons. The geometry points are con¬ 
nected via forces. In a first approximation, the elas¬ 
tic structures can be identified to represent the elastic 
cell membranes. However, more elastic structures can 
be added to the intra-/extracellular volume to mimic 
the visco-elastic properties of the cytoskeleton or the 
extracellular matrix. The user can implement biolog¬ 
ical mechanisms which operate on the cell represen¬ 
tations. For example, a new junction to a neighbor¬ 
ing cell might be created when the distance between 
two neighboring cell boundaries falls below a thresh¬ 
old distance. Similarly, a junction might be removed 
when overly stretched. 

6.2 Fluid and Fluid-Structure Inter¬ 
action 

The visco-elastic cell model represents the content of 
cells (the cytoplasm) as well as the extracellular space 
(the interstitial fluid and the extracellular matrix) 
as a viscous, Newtonian fluid. The intra- and ex¬ 
tracellular fluids interact with the elastic membrane, 
i.e. the fluids exert force on the membrane, and the 
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membrane exerts force on the fluids. Furthermore, 
the velocity field of the fluid, which is induced by 
the forces, moves and deforms the elastic structures. 
This interaction, well-known as fluid-structure inter¬ 
action, lies at the heart of the tissue model. Forces 
(e.g. membrane tension or cell-cell forces) acting on 
these points are exerted on the fluid by distributing 
the force to the surrounding fluid. Due to the local 
forcing, the fluid moves. At this step, the membrane 
point are advected passively by the fluid. As a result 
the forces need to be re-evaluated on the points. By 
repeating the forcing-advection steps, the interaction 
is realized iteratively. 

As a result of this iterative process, the (elastic) 
structures are coupled to the (viscous) fluid. De¬ 
pending on the parametrisation, this model allows 
to describe both elastic, or viscous, or visco-elastic 
material behaviour. The upper part of Figure [l] illus¬ 
trates the immersed boundary interaction. The im¬ 
plemented IB kernel function has bounded support, 
i.e. each geometry point influences and is influenced 
only its immediate neighborhood. Here, the dimen¬ 
sion of the kernel function is four by four (cf. Figure 
[l]). The fluid equations are solved using the Lattice 
Boltzmann method ([5]), which is described in de¬ 
tail in the supplementary material (section [To]) . The 
Reynolds number is typically 1, hence the regime 
is described by Stokes flovQ 

6.3 Signalling 

The signalling network is represented as a system 
of reaction-advection-diffusion processes. The elastic 
membranes may act as no-flux boundaries for com¬ 
pounds which only exist in the extra- or intracel¬ 
lular volume, respectively. The reaction-advection- 
diffusion solvers can be equipped with potentially 
coupled reaction terms in order to model signalling 
interactions of diffusing factors. Depending on the 
model, the signalling may impact the tissue dynam¬ 
ics. This can be done, for instance, by making the 

1 The Reynolds number reads Re = FT, with U be¬ 
ing a characteristic velocity, L a characteristic length scale, 
and v the kinematic viscosity. Assuming L = 10 —3 [m], 
U = 10 -8 [m/s] and v = 10 1 .. . 10 2 [m 2 /s], then Re = 

10 -13 .. . 10 -12 can be estimated l[T0)b 


mass source of the fluid dependent on the values of 
the reaction-advection-diffusion solvers such that the 
tissue expands locally (cf. Figure [l]). Furthermore, 
the diffusing compounds can be individually config¬ 
ured to diffuse freely across the entire domain, or only 
inside or outside the cells (e.g. using no-flux bound¬ 
ary conditions for the cell membranes). 

7 Software 

7.1 Cell Representation 

The cell geometries consist of two elements, the 
GeometryNodes, which act as the IB points, and the 
Connections, connecting pairs of GeometryNodes. 
A simplified cell is visualized in Figure [2j The 
Connections are attributed with a domainID flag, 
which is an identifier for the surrounded domain 
(respecting the counter-clockwise directionality con¬ 
vention). The domain identifier on the other side 
(on the right hand side) is zero by convention, 
representing the interstitial space. The domainID 
of the Connections are copied to the fluid and 
reaction-advection-diffusion solvers. Moreover, the 
domainID’s are associated with a cell type flag, 
cellType. By applying custom differentiation rules, 
the cellType of individual cells may be changed ac¬ 
cording to custom criteria; otherwise the all cells de¬ 
fault to cellType=l (with cellType=0 being the 
interstitial space, again). In this way, the reaction 
terms and the mass sources may be made dependent 
on specific cells, or specific cell types. 

7.2 User-provided Solvers 

The user can add the following routines: 
MassSolverXX, CDESolverXX, and BioSolverXX 
(XX being a name to be chosen). The MassSolverXX 
- as described above - adds or substracts mass 
from/to the fluid solver. The CDESolverXX is used 
to implement the reaction terms of the signalling 
models. Finally, the BioSolverXX can be used 
to execute biologically motivated operations on 
the geometry and the forces. Such an operation 
might be cell division, which is discussed in more 
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Figure 2: Elements of the Geometry Represen¬ 
tation. The cells are closed polygons, consisting of 
geometry nodes (discs in the top part) and connec¬ 
tions (shaded boxes in the top part) between each 
two geometry nodes. Each connection stores two 
references to its preceding and successive geometry 
nodes, and vice versa each geometry node stores two 
references to its preceding and successive connection 
(visualized by aggregation arrows in the top part). 
Directionality of the polygon is counter clockwise by 
convention. Each geometry node has a unique, im¬ 
mutable node ID attribute, which is allocated inter¬ 
nally upon creation of a new geometry node. Each 
connection features a domainID attribute, which de¬ 
notes the domain identifier of the domain on the 
left hand side. The domain identifier on the right 
hand side is by definition zero, representing the ex¬ 
tracellular space. Using the domainID of the connec¬ 
tions, the domainID of the lattice nodes is automat¬ 
ically set (lower part). Additionally, each domainID 
is associated with a cellType. The behaviour of 
the MassSolverXX, BioSolverXX and CDESolverXX 
can be made dependent on the domainID and/or 
cellType attributes by the user. 




Figure 3: Simplified UML diagram of impor¬ 
tant Classes. The classes which have to be provided 
by the user are shaded. XX refers to an arbitrary 
solver name. A The SimulationRunner controls the 
execution of the simulation. The GeometryHandler 
has a collection of PhysicalNodes, representing the 
lattice, a collection of BoundaryNodes wich are wo¬ 
ven into the lattice, and a Geometry object. The lat¬ 
ter contains the cell’s geometric information, namely 
the GeometryNodes and the Connections. The 
GeometryNodes and the Connections each have two 
references of the preceding and successive elements, 
as also explained in Figure [2j BioSolverXX ob¬ 
tains references from the GeometryHandler and the 
ForceSolver to alters states accordingly. Similarly, 
the MassSolverXX obtains a reference to the lattice 
and adds mass sources to the fluid. B To imple¬ 
ment new custom routines, the user must inherit from 
provided base classes (from BioBaseSolver for bi¬ 
ologically motivated routines, from BaseCDESolver 
for reaction-advection-diffusion processes, and from 
BaseMassSolver for mass modifying routines) 


detail in Subsection |7.5.4| Figure [3]4 summarizes 
the most important classes and their interactions. 
The classes which are subject to customization are 
shaded. In order to add a new customized routine 
(e.g. a mass modifying solver MassSolverXX, a 
reaction-advection-diffusion solver CDESolverXX, or 


a biologically motivated solver BioSolverXX), the 
user needs to inherit from their respective virtual 
base classes (cf. Figure [3^3). Figure [4] visualizes 
the routines, which are called iteratively by the 
SimulationRunner (cf. Figure [3jA) . 
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Figure 4: Iterative Processing in the Solver At 

initiation, the library loads the user-provided con¬ 
figuration files (containing global simulation param¬ 
eters, initial geometry, initial forces). During each 
iteration, the library’s class SimulationRunner (cf. 
Figure HF successively calls the physical routines 
(the Lattice Boltzmann method to solve the fluid and 
react ion-ad vection-diffusion processes, and the im¬ 
mersed boundary method to solve the fluid-structure 
interaction) and the biological routines (biologically 
motivated re-arrangement of the geometry, modifica¬ 
tions of the forces, etc.). The current configuration 
and optionally the entire solver states can be saved 
at a chosen frequency. 


7.3 Input and Output 

The communication to the user is achieved via the 
loading and dumping of configuration files. A gen¬ 
eral configuration file contains the global simulation 
parameters, such as the simulation time, the do¬ 
main size, the fluid viscosity, and the diffusion co¬ 
efficients for the reaction-advection-diffusion solvers. 
The geometry points and the corresponding geomet¬ 
rical connections are stored in a geometry file. A 
third file contains the forces, including forces between 
a pair of geometry points, freely defined forces or 
spatially anchored points. The fluid and reaction- 
advection-diffusion solver states may be written ei¬ 
ther to .txt files or in .vtk format and can be post- 
processed with third-party software (e.g. Matlab or 


ParaView). Optionally, the solver states can be saved 
in a loadable format to resume the simulation. 


7.4 Physical Processes 

7.4.1 Viscous and Elastic Behaviour 


The viscous behaviour is implemented using a repre¬ 
sentation of an incompressible fluid (solved using the 
Lattice Boltzmann method, cf. supplementary mate¬ 
rial), which converges to the Navier-Stokes equation 
in the hydrodynamic limit. The fluid is solved on a 
regular Cartesian and Eulerian grid. The membranes 
are represented by sets of points, which are connected 
to form closed polygons. A variety of forces may act 
on the membrane nodes, such as e.g. membrane ten¬ 
sions (cf. Subsection 7.4.3). The interaction between 
the fluid and the elastic structures is formulated using 
the Immersed Boundary method (cf. supplementary 
material). The membrane points move according to 
the local fluid velocity field in a Lagrangian manner. 


7.4.2 Reaction-Diffusion of Biochemical 
Compounds 

The biochemical signalling can be described by sets 
of coupled reaction-diffusion partial differential equa¬ 
tions. Similar to the fluid equations, these equa¬ 
tions are solved on a regular Cartesian and Eulerian 
grid (solved using the Lattice Boltzmann method, cf. 
supplementary material). The concentrations of the 
compounds can be accessed by other solvers, for ex¬ 
ample to make other processes such as cell division 
dependent on signalling factors. The cell boundaries 
can be chosen to be either invisible to the diffusing 
compounds, or to be no-flux boundaries. To account 
for advection, the fluid velocity field is directly trans¬ 
ferred from the fluid solver since the fluid and the 
reaction-diffusion lattices coincide spatially. The cou¬ 
pling of the solvers is visualized in Figure [1] 

7.4.3 Forces 

Forces are an integral part of the simulation environ¬ 
ment. A force is always connected to a membrane 
point. Any type of conservative force (which can be 
derived from a potential) can easily be implemented. 
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Currently, the following types of forces are imple¬ 
mented: 

• spring force between two geometrical nodes 

• spring force between a geometrical node and a 
spatial anchor point 

• free force acting on a geometrical node 

• horizontally or vertically sliding force (thus en¬ 
forcing only the y or x coordinate, respectively) 

• constant force between two geometrical nodes 


Application examples include constant forces be¬ 
tween two geometrical nodes that can be used to 
model constant membrane tension, which leads to the 
minimization of a cells perimeter (discussed in Sub¬ 
section 7.5.2). Moreover, a geometrical point can dy¬ 


namically explore its local neighbourhood and estab¬ 
lish a force to another geometrical point from another 
cell , thus mimicking cell-cell junctions (discussed in 
Subsection 7.5.3). 


7.5 Biological Processes 

The biological solvers (BioSolver) accommodate the 
functionalities that are related to biological processes. 
These processes may be mostly related to modifica¬ 
tions of the forces and the geometry. The BioSolver 
has full access to the compound concentrations. Fur¬ 
thermore, it is aware of the cells, whose geometries 
are stored individually. This enables the BioSolver to 
compute cell areas and averaged or integrated com¬ 
pound concentrations. Since all cells are individually 
tagged, cell behaviour can be made dependent on cell 
identity. Additionally to the cell identity, cells also 
carry a cell type tag, which can be changed depend¬ 
ing on run-time conditions. This latter functionality 
can be used to model cell differentiation. 

Consider a cell division event as an example. Here, 
a division plane has to be chosen. The choice of its 
position and direction is subject to the user’s model: 
the cell division plane might be set perpendicular to 
the cell’s axis of strongest elongation. Next, the cell 
has to be divided, which requires the removal of the 


corresponding geometrical connections, and the in¬ 
sertion of new geometrical nodes and connections to 
close the divided cells. 

Note that the concentration fields of the com¬ 
pounds, as well as the velocity- and pressure fields 
of the fluid solver are not directly altered in the bio¬ 
logical module. 

7.5.1 Control of Cell Area 

Depending on the biological model of the user, the 
cell area has to be controlled. By assuming that a cell 
might change its spatial extent in the third dimen¬ 
sion, the area might shrink or expand as a response 
to forces exerted by its neighbouring cells, which can 
effectively be modelled as an ’area elasticity’. In the 
limiting case, the cell resists external forces, main¬ 
tains its area and only reacts with changes of the 
hydrostatic pressure. In general, to control the area 
of cells, the reference area for each cell needs to be 
adapted. The reference area acts as a set point for 
a simple proportional controller, i.e. the local mass 
source Sk in the cell k is proportional to the area dif¬ 
ference between the current cell area (£), and the 
set point area A®: 


S k = a (A° k - (i)) (1) 

where cu is a proportional constant. More advanced 
control methods, such as e.g. proportional-integral 
control methods, can be realized easily. 

This approach of controlling the cell area can also 
be used to let cells grow or shrink in a controlled 
way, i.e. a cell differentiating into an hypertrophic 
cell type may grow in volume. Implementing this 
process would be as simple as setting the new target 
area as set point area. The area controller will bring 
the cell close to its new area. 

7.5.2 Membrane Tension 

The definition of forces acting between pairs of mem¬ 
brane points allows for simulating the cell’s mem¬ 
brane tension. By default, a constant contracting 
force Fi with magnitude (p m is applied to every pair 
of neighbouring membrane points. Hence the result¬ 
ing force on membrane point i is composed of a force 
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pointing to its preceding membrane point i — 1 , and a 
force pointing to its successive membrane point i + 1: 


*r-<p n 


1 *&i 


*^i+l 


\ x i -1 




|**h+l 




( 2 ) 


This approach can be interpreted as an actively re¬ 
modelled membrane: when stretched, new membrane 
is synthesized in order to not increase the membrane 
tension on longer time scales (hours). On the other 
hand, excessive membrane is degraded to abide the 
membrane tension. Therefore, the membrane tension 
minimizes the cell’s perimeter. Since the intracellu¬ 
lar fluid (and thus the cell area) is conserved in the 
absence of neighbouring cells and active mechanisms 
(c.f. Subsection |7.5.1 ), the cell assumes a circular 
shape. On short time scales (seconds), the passive 
(non-remodelled) elastic membranes can be modelled 
by using Hookean spring potentials. The membrane 
tension will then be proportional to deviation from 
the resting membrane perimeter. In both cases, the 
membrane is flexible (i.e. has no bending stiffness); if 
bending stiffness should be required by the user, this 
can be easily realised in a custom BioSolver. 

The implementation of membrane tension needs to 
consider the geometry remeshing. Whenever a new 
membrane point is inserted, it needs to get connected 
to its neighbours instantly, because the cell will be 
overly stretched in the absence of membrane tensions. 
A membrane point’s forces need to be removed upon 
its removal. Algorithmically, this is realized by re¬ 
moving and reconstructing all membrane forces at 
every time step. At this point, the magnitude of the 
membrane tension can be made dependent on sig¬ 
nalling factors. 

BioSolverMembraneTension is an example of 
a class managing the membrane tensions with 
immediate remodeling, and 

BioSolverHookeanMembraneTension implements 
simple Hookean springs. 


7.5.3 Cell Junctions 

A cell can create cell junctions to neigh¬ 
bouring cells. In the simplest case, 

each membrane point i uses the function 
getGeometryNodesWithinRadiusWithAvoidance- 


Closest to get the closest membrane point j of 
another cell, which is within a predefined cut-off 
radius l max , or zero if there is no such membrane 
point. Once a candidate membrane point fulfils 
the criteria, a new Hookean force Fi with a spring 
constant k) and resting length Zq is created: 

pcj _ f ^ \xj-Xi\ (I X 3 — X i \ — Jo) if \ x j ~ X i\ < 

% y 0 else 

(3) 

The cell junction forces are regularly (potentially 
not at every time step) deleted and renewed, where 
the frequency of cell junction renewal might reflect 
the cell junction synthesis rate. 

The function getGeometryNodesWithinRadius- 
WithAvoidance returns all membrane points of an¬ 
other cell, which are within a predefined cut-off ra¬ 
dius; the returned list might be empty. This opens 
up the possibility to introduce randomness by choos¬ 
ing the membrane point randomly from the candidate 
list. The probability to create a junction might de¬ 
pend on the junction length: the shorter, the higher 
the probability to form a new junction. Also the 
removal of membrane points might be randomized, 
and the probability made dependent on the junction 
length, i.e. overly stretched junctions are removed 
with higher probability. Even the membrane point 
whose junctions shall be updated might be chosen 
randomly. Again, the number of updated membrane 
nodes per time reflects the cell’s limited cell junction 
synthesis activity. 

The membrane points are internally stored in an 
fast neighbor list data structure, which is well suited 
for spatial range queries. BioSolverCellJunction is 
an example of a class responsible for cell junctions. 

7.5.4 Cell Division 

The cell division functionality requires several steps. 
First, criteria will have to be defined which cells shall 
be divided. Criteria might be maximal cell area, max¬ 
imal spatial expansion, or biochemical signals. Once 
a cell committed for division, the cell division plane 
will have to be chosen. Again, how to chose the plane 
is subject to biological modelling. A frequently used 
rule is to use a plane defined by a random direction 








vector and the center of mass of the cell. However, 
different rules can be readily implemented, such as 
random directions drawn from non-uniform proba¬ 
bility distributions (which, in turn, can be controlled 
e.g. by signalling factor gradients) or division planes 
perpendicular to the longest axis (22)- I n a next 
step, the two membrane segments are determined 
which intersect with the division plane; this is imple¬ 
mented in getTwoConnectionsRandomDirection or 
getTwoConnectionsLongestAxis. These two mem¬ 
brane segments are subsequently removed, and two 
new membrane segments across the cell are intro¬ 
duced, leading to a cut through the mother cell. Fi¬ 
nally, a new domain identity number has to be given 
to one of the daughter cells; the other daughter cell 
inherits the domain identity number from the mother 
cell. The new domain identity number is set to the 
largest domain identity number plus one, and it is au¬ 
tomatically copied to the physical grid. Both daugh¬ 
ter cells by default inherit the cell type flag from the 
mother cell, which is also automatically copied to the 
physical grid. 

The basic cell division functionality is implemented 
in the class BioSolverCellDivision. 

7.5.5 Differentiation 

Differentiation changes the cell type flag of the 
cells according to user-defined, biologically motivated 
rules. These rules might be based on the cell area, 
or on a signalling factor concentration, possibly in¬ 
tegrated over the cell area. Once being committed 
for differentiation, the cell changes its cell type flag 
according to the rule. The new cell type flag will 
be automatically copied to the physical grid. The 
cell type flag can be used to make signalling dynam¬ 
ics, but also other biologically motivated processes 
dependent on the cell type. 

The association between the domain identifier 
flags and the cell type flags is stored in the 
cellTypeTrackerMap_, which is a member of the 
GeometryHandler. This makes sure that all 
BioSolverXX classes have easy access to this infor¬ 
mation. A basic implementation of the differentiation 
control can be found in BioSolverDif f erentiation. 


7.6 Accuracy and Performance 

The Lattice Boltzmann schemes are second or¬ 
der accurate, and the explicit immersed boundary 
method is first order accurate in space and time. 
The internal data structure uses a fast neighbor 
list (cell list) implementation to optimize for range 
queries (e.g. searching for other cells in the lo¬ 
cal neighborhood), which exhibits a search complex¬ 
ity of O (N ), with N being the number of mem¬ 
brane points to represent the cells. Many iter¬ 
ative computations (LB and IB routines such as 
particle streaming and collision, gathering of ve¬ 
locity and scattering of force) are parallelized us¬ 
ing the shared memory paradigm. However, a few 
computational steps cannot be parallelized. This 
is typically the case when write-operations occur 
on shared data structures, such as the data struc¬ 
tures storing the geometry nodes and the force 
structs (e.g. in ForceSolver: :deleteForceType() 
and GeometryHandler:: compute Areas ()). More¬ 
over, the geometry remeshing (refining and coarsen¬ 
ing) functions as well as the data I/O are not parel- 
lelized, but are assumed to occur much less frequently 
than the actual fluid and reaction-advection-diffusion 
solvers. Therefore, since the fraction of sequential 
code is not negligible, the software should best be 
run on fast multi-core processors. 


7.7 Tools, Dependencies and Docu¬ 
mentation 

A compiler with C++0x support (such as GCC 

4.7 or higher) is required. The software de¬ 
pends on Boost (http://www.boost.org; 1.54.0 or 
higher), OpenMP, CMake (http://www.cmake.org) 
and vtk (http://www.vtk.org/; 5.8 or higher). The 
source code is extensively documented using Doxygen 
(http://www.stack.nl/~dimitri/doxygen). Git 
(http://git-scm.com) is used for version control. 
The software has only been tested on linux operating 
systems. 
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7.8 Availability 

The documentation and source code are available on 

http://tanakas.bitbucket.org/lbibcell/ 
index.html 


ken if overly stretched, and new junctions are formed 
(according to Eq. © ). At the boundary of the tis¬ 
sue, the cells try to reach a spherical shape, while in 
the middle mainly characteristic penta- and hexago¬ 
nal shapes emerge (cf. Figure [5^ and Supp. 10.4). 


8 Application Examples 

8.1 Cell Division, Differentiation and 
Signalling 

To demonstrate the capabilities of the software we 
first consider a tissue model with cell-type specific 
cell division and signalling-dependent differentiation 
(Figure [5]). In the beginning, a circular cell with ra¬ 
dius R = 10 is placed in the middle of a quadratic 400 
by 400 domain (Figure [ 5 }^). Iso-pressure boundary 
conditions are set at the border of the domain. The 
initial cell is of red cell type, which is proliferating at 
a high rate. When considering a single layer epithe¬ 
lium, mass uptake, which is needed for modelling cell 
growth and finally proliferation, is assumed to occur 
from the apical cavity through the apical membrane. 
Additionally, the initial cell secretes a signalling fac¬ 
tor X which inhibits differentiation of the red cell type 
into the green cell type. Once the cell area doubled, 
the cell is divided in a random direction (cf. Figure 
[5j3). The daughter cells inherit the cell type, but only 
the mother cell continues to express the signalling 
molecule X. All cells of red type integrate the con¬ 
centration of X over their area. For low signalling 
levels, the red cell type differentiates into the green 
cell type. The green cell type does not grow and only 
divides if external forces stretch the cell. In Figure 
Up, the daughter cell’s signalling level dropped after 
cell division, and differentiation occurred. After sev¬ 
eral rounds of cell division, a tissue starts to form (cf. 
Figure|5jD). The cells close to the secreting initial cell 
remain protected from differentiation, whereas more 
distant cells differentiate irreversibly. Due to the ran¬ 
domly chosen cell division axis, it might happen that 
the proliferating red cells get trapped (cf. Figure [5^)). 
The expression of X is switched off at time t=5000, 
thus leading to complete differentiation shortly after 
(cf. Figure^). After proliferation stopped, the cells 
slowly rearrange because cell-cell junctions are bro- 


8.2 Turing Patterning on Growing 
Cellular Domains 

To demonstrate the importance to investigate mor- 
phogenic signalling hypotheses on dynamically grow¬ 
ing domains with cellular resolution, we solved a 
reaction-diffusion system, featuring the well-known 
diffusion-driven Turing instability ([4l]), on a prolif¬ 
erating tissue. Figure [ 6 ]A illustrates the interaction 
between a ligand L and its receptor R. Here, we as¬ 
sume that one ligand dimer molecule L binds to two 
receptors R, forming the complex R 2 L which induces 
upregulation of the receptor on the membrane (e.g. 
0 ). Unbound receptor is turned over at a linear rate. 
The ligand can diffuse freely across the tissue and the 
entire domain, whereas the diffusion of the receptor 
is limited to a single cell’s apical surface and is much 
slower. The dynamics can be formulated as a system 
of non-dimensional partial differential equations: 

d t R = AR + 7 (a - R + R 2 L) (4) 

d t L = dAL + 7 (b - R 2 L) (5) 

where 7 is a reactivity constant, a and b produc¬ 
tion constants, and d the relative diffusion coefficient 
of ligand and receptor. We note that the equations 
correspond to the classical Schnakenberg-type Turing 
mechanism mmi It has previously been shown 
that such a receptor-ligand interaction can explain 
symmetry breaking in various morphogenetic systems 

(eh m m [201 sni mi). 

Depending on the type of domain we observe very 
different patterns. On a continuous domain we obtain 
the well known regular spot pattern (Figure [6^3). On 
an idealized static cellular domain an overall regular 
pattern with irregular internal structure (Figure [ 6 p) 
can be observed. Decreasing the simulation parame¬ 
ter 7 , which inversely controls the distance between 
the spots, leads to even more unexpected patterns: 
for 7 = 100 , the local regularity is completely lost 
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(Figure [ 6 jD ). Finally, on a dynamically growing cel¬ 
lular domain, where the local proliferation rate was 
set proportional to the R 2 L signal, we obtain irregu¬ 
lar patterns (Figure | 6 jE). For a lower value 7 = 100, 
clusters of cells with high R 2 L signalling levels emerge 
(Figure | 6 ^). In conclusion, even relatively simple sig¬ 
nalling mechanisms can lead to significantly different 
results, depending on how the tissue is represented. 

9 Discussion 

We developed an extendible and open-source cell- 
based simulation environment, which is tailored to 
study morphogenetic problems. The novel frame¬ 
work permits the coupled simulation of a physically 
motivated visco-elastic cell model with regulatory sig¬ 
nalling models. Processes such as viscous dissipation, 
elasticity, advection, diffusion, local reactions, local 
mass sources and sinks, cell division and cell differ¬ 
entiation are implemented. By applying our frame¬ 
work to Turing signalling systems we show that the 
signalling systems may behave very differently on dy¬ 
namic tissues than on simple continuous tissue rep¬ 
resentations. We therefore advocate to test continu¬ 
ous morphogenetic signalling models on dynamically 
growing cellular domains. 

The presented framework permits to study a vari¬ 
ety of mechano-regulatory mechanisms. By making 
the cell division orientation dependent on signalling 
cues, the effect on the macroscopic tissue geometry 
may be studied. Cell migration can be modelled by 
introducing gradient-dependent forces on specific cell 
types. Cell sorting may be achieved by specifying 
multiple cell types with differential cell-cell junction 
strengths. The framework is specifically designed to 
study the mutual effects of signalling and biophysical 
cell properties. 

The visco-elastic cell model represents cell shapes 
at very high resolution and is thus, unlike the ver¬ 
tex model, not restricted to densely packed tissues. 
Furthermore, hydrodynamic interaction, membrane 
tension and hydrostatic pressure are integral compo¬ 
nents of the model. The fact that a velocity field 
is available on the entire domain is a critical advan¬ 
tage to account for advection of the signalling com¬ 


ponents, thus allowing for a spatial description of in¬ 
tracellular concentrations. The model is, however, 
not easily extendable to the third dimension. Since 
a meshing of the surface will be required, the algo¬ 
rithmic and computational complexity are expected 
to be significant and subject to future work. The 
presented framework is, however, ideal to study in¬ 
trinsically two-dimensional morphogenetic problems, 
such as apical surface dynamics of epithelia as studied 
previously also by [ 8 ] and [17! in 2D. 
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Figure 5: Cell Division, Differentiation and Sig¬ 
nalling. A The initial configuration consists of a sin¬ 
gle, circular cell of type red. The red cell type pro¬ 
liferates at a high rate. The initial cell is tagged and 
expresses a signalling molecule X which inhibits dif¬ 
ferentiation. B The first cell division occurs. The di¬ 
vision axis is chosen randomly. The daughter cell in¬ 
herits the cell type from the mother cell, but only the 
mother cell keeps expressing the signalling molecule 
X. C The signalling level (the spatially integrated 
concentration of the signalling molecule) drops in 
cells far away from the initial cell and differentia¬ 
tion into the green cell type occurs. The green cell 
type does not grow intrinsically, and only divides if 
overly stretched by external forces. D The highly 
proliferating red cells are trapped in the forming tis¬ 
sue due to the randomly chosen cell division axis. 
At t = 5000, the expression of the differentiation in¬ 
hibiting molecule X is switched off, which leads to 
the differentiation of the remaining red cells. E In 
the absence of high proliferation, the cells rearrange 
to maximize the perimeter/area ratio. Characteristic 
penta- and hexagonal cell shapes emerge (cf. Supp. 
10.4). Cells close to the boundary try to take a cir¬ 


cular shape. 
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Figure 6 : Turing Patterning on Growing Cellu¬ 
lar Domains A Turing instability can be achieved 
by Schnakenberg-type reactions, involving a slowly 
diffusing compound R, here interpreted as a receptor, 
and a fast diffusing compound L, here interpreted as 
a freely diffusing ligand. One ligand molecule binds 
to two receptors, leading to the complex R 2 L. The 
complex can be interpreted as a biological signal. B 
The model is solved on a continuous square lattice 
(using d = 1, 7 = 800, a = 0.1, b = 0.9), result¬ 
ing in the classical regular spot-pattern. The bio¬ 
logical signal R 2 L is shown. C The same system as 
in B is solved on an idealized static cellular domain, 
i.e. the diffusion of the receptor R is restricted to a 
cell. The emerging biological signal R 2 L is now dis¬ 
tributed irregularly. D The same system as in C, but 
with 7 = 100 , is solved on a idealized static cellular 
domain. Fewer cells show significant levels of signal 
R 2 L and no regular pattern can be found (salt-and- 
pepper pattern). E The same system as in C is solved 
on a growing cellular domain. The proliferation rate 
of a cell is set proportional to its signal R 2 L. The re¬ 
sulting pattern features regularity on a larger scale, 
but the local patterning significantly differs from the 
behaviour on continuous (B) and static cellular (C) 
domains. F The same system as in D is solved on 
growing cellular domain. The proliferation rate of a 
cell is set proportional to the local intensity of the 
signal R 2 L. Clusters of active cells with high levels 
of R 2 L emerge. 
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10 Supplementary Material 

10.1 Lattice Boltzmann Method 

For the fluid, the standard Lattice Boltzmann scheme 
(with the single-relaxation-time Bhatnagar-Gross- 
Krook collision operator) is used ([5]). It has been 
shown, that the incompressible Navier-Stokes equa¬ 
tions can be recovered in the hydrodynamic limit 
(0). The Boltzmann equation is discretized on a 
D 2 Q 10 lattice and reads: 


fi(x + Vi,t + 1 ) - fi(x,t) = 

- - (.fi (x, t) - ff q (x, t)) + A t^Y- ( v i ■ /) ( 6 ) 

r cf 

where fi denotes the particle distribution function in 
the direction i, and can be interpreted as the proba¬ 
bility density of finding particles with velocity Vi at 
time t at position x. The relaxation time r is re¬ 
lated to the kinematic viscosity v = c 2 (r — |). For 
isothermal flows, the speed of sound c s is defined 
as c s = c/y/3 using the lattice speed c = Ax/At. 
The lattice spacing and the time step are chosen as 
Ax = At = 1 to guarantee consistency with the lat¬ 
tice. The last term of Equation ^ represents the 
external body force am)- 

Equation © implies a two step algorithm: firstly, 
the distribution functions fi perform a free flight to 
the next lattice point (left hand side). Secondly, on 
each lattice point, the incoming distribution func¬ 
tions collide and relax towards a local equilibrium 
distribution f/ q , which is controlled by the relaxation 
time r (right hand side). The equilibrium distribu¬ 
tion is taken as: 


/f = u iP 


Vi • u 


(ViVi — u 2 ) : uu 

2 ^ 


(7) 


For the D2Q9 lattice, the population weights uji can 
be found as ujq = 4/9, Ccq _4 = 1/9, and = 1/36. 

Finally, the macroscopic quantities (density p and 
momentum density pu) can be computed using the 
zeroth and first order moments: 


p = Ef> 


8 

p u = J2 f iVi 

i =0 


The fluid pressure p is related to the mass density p 
as p = pc 2 . 

For solving the reaction-advection-diffusion equa¬ 
tion of a compound 0 , a multi-distribution function 
approach is chosen mm ), i.e. an additional Lattice 
Boltzmann solver is coupled to the fluid solver. We 
implemented different schemes (D2Q4,D2Q5), which 
exhibit slightly different stability and accuracy am). 
For the D2Q5 [^] scheme, we follow m and m and 
write for the Lattice Boltzmann equation: 

gi(x + Vi,t + 1 )-gi ( x , t) = —— (gi (x, t ) - g\ q (x, t)) 

td 

(9) 

The equilibrium distribution functions is taken as: 


gT = <M 


1 + ~2 (Vi ■ u) 


( 10 ) 


Instead of the local first moment, the velocity field 
u from the fluid is transferred. The weights Wi are 
chosen as wo = 1/3 and 4 = 1/6. The relaxation 
parameter td is related to the diffusion coefficient as 
D = c 2 (td — |), and the local compound density 0 
reads: 

4 

i =0 


with the fluid velocity u and the fluid density p. The 
operator : denotes the dyadic tensor scalar product. 

2 the nine-velocity lattice in two dimen¬ 
sions is defined as vq = [0,0], v% = 

[cos (7r (i — 1)/2), sin (7r (i — 1)/2)] for i = (1,2, 3, 4}, 

and Vi = [cos (tt (i — 1) /2) + 7r/4, sin (n (i — 1) /2 + 7r/4)] 
for i = (5, 6, 7, 8} 


All variables are expressed in Lattice Boltzmann 
units Sx,St and have to be converted to physical 
units. 


3 the five-velocity lattice in two dimensions is defined as 
vo = [0,0] and Vi = [cos (tt (i — 1) /2), sin (tt (i — 1) /2)] for 
* = { 1 , 2 , 3 , 4 } 
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10.1.1 No-flux Boundary Condition for 
Reaction-Advection-Diffusion Equa¬ 
tions 

The missing incoming distribution functions are ap¬ 
proximated by equilibrium distributions. The mo¬ 
mentum is spatially first order interpolated between 
the fluid in direction of the missing population, and 
the known zero momentum at the wall. The density 
is spatially first order extrapolated from the fluid. 


10.1.2 Pressure Boundary Condition for the 
Fluid Equations 

Within the computational domain, lattice points can 
act as internal pressure boundaries, i.e. points with 
prescribed fluid pressure, which may be needed to 
implement mass sinks in case of growing cells/tissues. 
All distribution functions are rescaled such that the 
prescribed pressure is obtained. The velocity field is 
not affected, since the rescaling factor 7 cancels: 


ecuted by integrating over the material coordinates: 
p(x,t) = j M (g, r, s,t)S(x-X (g, r, s, t)) dqdrds (13) 

f (x,t) = J F (g, r, s,t)S(x-X (g, r, s, t)) dqdrds (14) 

where S (•) denotes the delta Dirac function. For the 
numerical implementation, the Delta Dirac function 
is approximated in such a way that it covers multi¬ 
ple grid points, on which the conversions in Equation 
(13) and (15) are evaluated. The equation of motion 
of the Lagrangian particles reads: 

dX , 


u 


(X (q, r,s,t) ,t) as- J u (x, t) 5 (x - X ( q, r, s, t)) 


dx 


(15) 


Using the Hodge decompositioiQ the material i 
described by 


is 


u (x, t) 


Tn 7 

tiiUM 


( 12 ) 


10.2 Immersed Boundary Method 

The PDE’s are prefentially solved on a Eulerian 
grid, whereas cells are naturally represented in a La¬ 
grangian manner. The two frames can be bridged 
using the immersed boundary method, which was in¬ 
troduced in [28] . This technique is appropriate for 
complex and moving boundaries and fluid-structure 
interactions. The material equations are solved on 
a Eulerian grid, whereas the boundary equations are 
expressed in a Lagrangian way. An introduction can 
be found in [24] . and detailed discussions in [29.. 

The curvlinear (boundary) coordinates (g, r, s) are 
attached to a material point. At time £, the coor¬ 
dinates in the Eulerian framework read X (g, r, s, t). 
The material derivative ^ (sc, t) = (g, r, s, t) is 

the acceleration of whatever material point is at po¬ 
sition x at time t. The conversion from Lagrangian 
variables (e.g. mass density M (g, r, s, t) or force den¬ 
sity F (g, r, 8 , t)) to their Eulerian counterparts is ex- 


P 


du 

dt 


+ (u • V) u 


V • u 


—Vp + pA u + / (16) 
0 (17) 


To couple the Immersed Boundary method to the 
Lattice Boltzmann method, the LB forcing term f 
has to be obtained from the Lagrangian force F fol¬ 
lowing a similar way as [9]. The Lagrangian force 
Fi of the geometry point i is composed of membrane 
forces (cf. Equation ©) and cell junction forces (cf. 
Equation ([2])): 

Fi = F™ + F? (18) 


Fi is distributed to the local fluid neighborhood ac¬ 
cording to Equation ( [14] ) , and the Eulerian force f 
on the lattice point x is then used in the LB collision 
term in Equation ©• 

The discretized delta Dirac function for one spatial 
dimension is defined as: 



if ||r|| < 2 
if |H| > 2 


(19) 


4 any vector field u can be written as — f = —Vp + w, 
where V • w = 0 
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and the two-dimensional discretized delta Dirac func¬ 
tion is a multiplication of Equation 

S (r ) = S 1D ( r ) S 1D ( r ) (20) 

This implies that the Lagrangian force of a mem¬ 
brane point i is distributed to the 4 x 4 = 16 nearest 
fluid lattice points. By discretizing the membrane 
into mebrane points i, the discretized delta Dirac 
functions overlap. 

10.3 Validations 

To validate the implementation of the fluid solver, the 
reaction-advection-diffusion solver and the immersed 
boundary method, validation tests using problems 
with known analytical solutions are executed. 

The dimensionless problem is as follows: a circular 
cell (diameter 1) is placed in the middle of a square 
domain of size 10 by 10. The end simulation time is 
T = 1. Inside the cell, a diffusing species is defined, 
where the cell membrane acts as a no-flux bound¬ 
ary condition. No production and degradation of the 
species occurs. The initial condition of the species 
is set to be uniformly 1. The diffusion coefficient as 
well as the kinematic viscosity of the fluid are set to 
1. No membrane forces are applied; iso-pressure out¬ 
flow boundary conditions are applied to the domain 
boundaries. Please see movieS3.avi for an impres¬ 
sion of how the concentration decreases as the area 
increases. 

Three different scenarios are considered: 

• Case I: a constant point mass source is placed in 
the middle of the cell (dimensionless strength is 
1); see Fig. [7^ 

• Case II: a constant uniform mass source inside 
the cell (dimensionless strength is log( 2 )/A 0 , 
where = 0 . 5 2 7 r is the initial cell area; this 
leads to a doubling of the cell area in time T), 
see Fig. [7j3 

• Case III: a uniform mass source inside the cell, 
where the strength is scaled linear proportional 
to the species concentration (= 1 • C(t), where 
C(t) is the concentration of the species; this 



Figure 7: Validation Cases: Evolution of Area 
and Concentration T is the non-dimensional nor¬ 
malized time, A the non-dimensional normalized 
area, C the non-dimensionalized normalized concen¬ 
tration. The column XI shows area A vs. time T; 
column X2 shows concentration C vs. time T; and 
column X3 shows concentration C vs. area A. A 
Case I: constant point-like mass source; B Case II: 
constant uniform mass source; C Case III: uniform 
mass source proportional to the concentration C. 


leads to a linear mass/area increase in time), see 

Fig. [7p 

For all three cases, the circular cell starts to grow, 
and the species is diluted, i.e. the species concentra¬ 
tion drops as the area increased. Ideally, the inte¬ 
grated species concentration is conserved. The rel¬ 
ative errors of the simulations w.r.t. the analytical 
solutions are computed for both the area and con¬ 
centration evolution. 

The simulations are executed at different levels of 
space and time discretization. As a measure, we take 
the diameter resolution of the circular cell, which is 
set to { 10 , 20,40, 80}. As a result, the discretized val¬ 
ues (LB denotes Lattice Boltzmann units, and ND 
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denotes non-dimensional) for all cases can be found 
in Table Q] 


circle diameter 

10 

20 

40 

80 

domain size 

100x100 

200x200 

400x400 

800x800 

t fluid 

1 

1 

1 

1 

v 1 " 13 fluid 

1/6 

1/6 

1/6 

1/6 

r diff 

1 

1 

1 

1 

V 713 diff 

1/6 

1/6 

1/6 

1/6 


0.1 

0.05 

0.025 

0.0125 

rj~iL/13 

le4 

4e4 

1.6e5 

6.4e5 

S t 

le-4 

25e-6 

625e-8 

15625e-10 

s LB =S ND •% 

-X 

le-2 

le-2 

le-2 

le-2 


Table 1: Lattice-Boltzmann discretization of Cases I- 
III. ND and LB denote non-dimensional and Lattice 
Boltzmann units, respectively 

The validation results for Case I are given in Fig¬ 
ure [8j Figure [8]A shows the relative error tempo¬ 
ral evolution of the cell area for the different resolu¬ 
tions, and [8^3 the relative error temporal evolution 
of the integrated concentration. The concentration 
error is much more noisy because it is only evaluated 
on lattice points. When the cell grows, the mem¬ 
brane sweeps across the lattice, and whenever a new 
lattice site enters the intracellular domain, the inte¬ 
grated concentration increases abruptly. To compute 
the convergence w.r.t. to spatial and temporal resolu¬ 
tion, the relative errors at T ND = 1 are evaluated for 
the area, but the temporal mean is taken for the con¬ 
centration convergence analysis. Figure[8]C shows the 
convergence analysis for the area. For the finer grids, 
the convergence order is approximately 1. Figure [8p 
depicts the convergence analysis for the integrated 
concentration, where the fitted order of convergence 
is « 1.2. 

The validation results for Case II are given in Fig¬ 
ure [9j Figure shows the relative error temporal 
evolution of the cell area for the different resolutions, 
and [9}3 the relative error temporal evolution of the 
integrated concentration. The concentration error is 
much more noisy because the it is only evaluated on 
lattice points. When the cell grows, the membrane 
sweeps across the lattice, and whenever a new lattice 
site enters the intracellular domain, the integrated 
concentration increases abruptly. To compute the 
convergence w.r.t. to spatial and temporal resolu¬ 
tion, the relative errors at T ND = 1 are evaluated for 


-DIA=10 

- DIA=20 

u.uo 
-D 0.045 


-DIA=10 

- DIA=20 

- DIA=40 

- DIA=80 

0.04 

: -0.035 


- DIA=40 

- DIA=80 


0.2 0.4 0.6 0.8 

t/T 



0.2 0.4 0.6 0.8 

t/T 



D 



log (resolution) 


log (resolution) 


Figure 8: Case I Validation: Point Source. An 

initially circular cell is growing through addition of 
mass. The mass is added by using a point source 
in the middle of the domain. Four different spatial 
resolutions have been computed, corresponding to a 
resolution of DIA = {10, 20,40, 80} lattice points for 
the diameter of the initial circular cell. A The time 
evolution of the relative error of the cell area for 4 
different lattice resolutions. B The time evolution 
of the relative error of the spatially integrated con¬ 
centration (mass conservation) for 4 different lattice 
resolutions. C Convergence plot of the relative error 
of the cell area as a function of the lattice resolution. 
The fitted slope is -1.80. The dashed line denotes a 
slope of -2. D Convergence plot of the relative error of 
the integrated concentration (mass conservation) as 
a function of the lattice resolution. The fitted slope 
is -1.22. The dashed line denotes a slope of -1. 


the area, but the temporal mean is taken for the con¬ 
centration convergence analysis. Figure [9]C shows the 
convergence analysis for the area. For the finer grids, 
the convergence order is approximately 1. Figure [9^3 
depicts the convergence analysis for the integrated 
concentration, where the fitted order of convergence 
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log (resolution) 


Figure 9: Case II Validation: Uniformly Dis¬ 
tributed Source. An initially circular cell is grow¬ 
ing through addition of mass. The mass is added by 
applying a constant uniform mass source on the en¬ 
tire cellular domain. Four different spatial resolutions 
have been computed, corresponding to a resolution 
of DIA = {10, 20,40, 80} lattice points for the diam¬ 
eter of the initial circular cell. A The time evolution 
of the relative error of the cell area for 4 different 
lattice resolutions. B The time evolution of the rel¬ 
ative error of the spatially integrated concentration 
(mass conservation) for 4 different lattice resolutions. 
C Convergence plot of the relative error of the cell 
area as a function of the lattice resolution. The fit¬ 
ted slope is -0.99. The dashed line denotes a slope 
of - 1 . D Convergence plot of the relative error of 
the integrated concentration (mass conservation) as 
a function of the lattice resolution. The fitted slope 
is -1.43. The dashed lines denote slopes of -1 and -2, 
respectively. 


is « 1.4. 

The validation results for Case III are given in Fig¬ 
ure [TO} Figure EUK shows the relative error tempo¬ 
ral evolution of the cell area for the different resolu¬ 
tions, and [Top the relative error temporal evolution 


of the integrated concentration. The concentration 
error is much more noisy because the it is only eval¬ 
uated on lattice points. When the cell grows, the 
membrane sweeps across the lattice, and whenever a 
new lattice site enters the intracellular domain, the 
integrated concentration increases abruptly. To com¬ 
pute the convergence w.r.t. to spatial and temporal 
resolution, the relative errors at T ND = 1 are evalu¬ 
ated for the area, but the temporal mean is taken for 
the concentration convergence analysis. Figure Hop 
shows the convergence analysis for the area. For the 
finer grids, the convergence order is approximately 
1. Figure [Top depicts the convergence analysis for 
the integrated concentration, where the fitted order 
of convergence is ~ 1.5. 

Finally, the immersed boundary force distribution 
is validated (Case IV). The initial setup is similar 
to the Cases I-III, but a constant membrane force 
is applied, and there is no mass source. Due to nu¬ 
merical leakage, the cell starts to shrink and looses 
mass/area. The dimensionless force is chosen to be 
F nd = le3, and the dimension conversion is de¬ 
scribed in Table [2j The discretization of the vali¬ 
dation setup can be found in Table [3] 



ND 

LB 

mass density 

p ND 

p lb 

mass units [M] 

[i] 

[ 

L> j 

mass density units 2D [M / L 2 ] 

[i] 


1 rLB-t-N D 

m 

force 

p N1J 

F LB 

force units [ ML/T 2 ] 

[i] 


cLB-t-NDcLB—t-ND 

m °x 

(gLB^ND^ J 


Table 2: Lattice-Boltzmann discretization of mass 
and force. ND and LB denote non-dimensional and 
Lattice Boltzmann units, respectively. 

Besides the spatial and temporal discretization of 
the lattice, also the discretization of the immersed 
boundary is studied in Case IV. The distance be¬ 
tween any two neighboring boundary points is shorter 
than the parameter MAXLENGTH, which is always a fac¬ 
tor of the lattice discretization S x . If MAXLENGTH is 
set to 0.5, then the distance between two boundary 
points does not exceed half of the lattice spacing. 
MAXLENGTH=0.5 is a frequently used constant ([3D]); 
however, here, we also tested smaller (0.1) and larger 
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Figure 10: Case III Validation: Fully Coupled 
System. An initially circular cell is growing through 
addition of mass. The mass is added by applying a 
mass source on the entire cellular domain, but the 
mass source strength is proportional to the concentra¬ 
tion of the intracellular concentration. Therefore, the 
cellular area should increase linearly in time. Four 
different spatial resolutions have been computed, cor¬ 
responding to a resolution of DIA = {10,20,40,80} 
lattice points for the diameter of the initial circular 
cell. A The time evolution of the relative error of 
the cell area for 4 different lattice resolutions. B The 
time evolution of the relative error of the spatially 
integrated concentrations (mass conservation) for 4 
different lattice resolutions. C Convergence plot of 
the relative error of the cell area as a function of the 
lattice resolution. The fitted lattices is -1.01. The 
dashed line denotes a slope of - 1 . D Convergence 
plot of the relative error of the integrated concentra¬ 
tion (mass conservation) as a function of the lattice 
resolution. The fitted slope is 1.50. The dashed lines 
denote slopes of -1 and -2, respectively. 


(2.5) values. 

The validation results of Case IV are showed in 


circle diameter 

80 

160 

320 

640 

domain size 

100x100 

200x200 

400x400 

800x800 

t fluid 

1 

1 

1 

1 

v 1 ' 1 * fluid 

1/6 

1/6 

1/6 

1/6 

t diff 

1 

1 

1 

1 

V 73 diff 

1/6 

1/6 

1/6 

1/6 

fix 

0.1 

0.05 

0.025 

0.0125 

rj-iLi fci 

le4 

4e4 

1.6e5 

6.4e5 

fit 

le-4 

2.5e-5 

6.25e-6 

1.5625e-06 

mass density p Nu = 1 

1 

1 

1 

1 

“7 p 1 ^ /c X2 

fim ~ (<**) 

le-2 

2.5e-3 

6.25e-4 

1.5625e-4 

piLB pND 

_ s m s x _ 

le-3 

5e-4 

2.5e-4 

1.25e-4 


Table 3: Lattice-Boltzmann discretization of the con¬ 
vergence analysis of Case IV. ND and LB denote 
non-dimensional and Lattice Boltzmann units, re¬ 
spectively. 


Figure [TT| In Figure [lT]A, the evolution of the rela¬ 
tive errors of the area are shown for different lattice 
resolutions. Figure \n£ shows the convergence as a 
function of MAXLENGTH. The MAXLENGTH= {0.1, 0.5} 
values lead to almost similar results, thus confirming 
that MAXLENGTH=0.5 is indeed sufficient. The conver¬ 
gence as a function of the lattice resolution is shown 
in Figure [TT|C. The fitted order of convergence is ap¬ 
proximately 1. MAXLENGTH= {0.1, 0.5} lead to indis¬ 
tinguishable results; however, even MAXLENGTH=2.5 
still leads to converging results. 

Here, for the convergence analyses, the fluid, mass 
sources, fluid outlet boundary conditions, the advec- 
tion and diffusion of a species, its no-flux boundary 
condition, the immersed boundary and membrane 
forces are considered. We showed that all aspects are 
converging and conclude that the net order of conver¬ 
gence of the fully coupled system is 1. Depending on 
the aspect under consideration, it can be higher. It 
is well known that the standard Immersed Boundary 
method suffers from fluid leakage, because the dis¬ 
cretized kernel function is not divergence-free. Sev¬ 
eral improvements have been proposed (e.g. [[301126] ). 

10.4 Comparison to Farhadifar et al., 
2007 

The cell topologies are compared to results of the 
cell vertex model as presented in [8]. Starting from a 
single cell, a tissue with more than 1000 cells is grown 
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Figure 11: Case IV Validation: Mass Conser¬ 
vation A constant membrane tension is applied to 
an initially circular cell. The contractile force will 
lead to a shrinking cell due to mass leaking of the 
immersed boundary. Here, the influence of the im¬ 
mersed boundary discretization is investigated. Four 
different spatial resolutions have been computed, cor¬ 
responding to a resolution of DIA = {10,20,40,80} 
lattice points for the diameter of the initial circu¬ 
lar cell. MAXLENGTH denotes the maximally toler¬ 
ated distance between each two immersed boundary 
points (measured in space discretization units). A 
The time evolution of the relative error of the cell 
area for 4 different lattice resolutions, and for 3 val¬ 
ues of MAXLENGTH. B The relative error of the cell 
area as a function of MAXLENGTH for different lattice 
resolutions. C Convergence plot of the relative error 
of the cell area as a function of the lattice resolution, 
plotted for 3 values of MAXLENGTH. The fitted slope of 
the two finest lattices is 0.9891. 


an area of 380 (in LB units). The membrane tension 
is varied ( F LB = {0.001,0.01,0.02} in LB units). 

The snapshots in Figure [T2| show the topology for 
different values of membrane tension (A very low 
membrane tension; B medium membrane tension; 
C high membrane tension). The higher the mem¬ 
brane tension, the higher the rearrangement activity 
because cell-cell junctions are broken more easily, and 
the cells can more easily assume a preferential shape 
(hexagonal in the domain, circular at the boundary). 
In Figure [T2p, the tissue shown in figure & is re¬ 
laxed, i.e. cell division is stopped, and only rear¬ 
rangement is occuring towards an equilibrium config¬ 
uration. 

For these cases, the number of vertices are counted. 
The result is shown in Figure [l3j together with the 
experimental and model results of [8]. Although the 
topologies for different membrane tensions differ sig¬ 
nificantly, the distribution of the number of vertices 
(equivalent to the number of neighbors) is not sig¬ 
nificantly affected. The relaxed topology, however, is 
shifted towards higher occurrence of hexagons. Over¬ 
all, our results are in agreement with the experimen¬ 
tal and simulation results obtained by [A. 


with a constant mass source of 0.004 (in LB units). 
A cell is divided in a random direction if it exceeds 
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Figure 12: Cell Topologies for Different Mem¬ 
brane Forces. A The membrane force constant 
is small ( F LB = 0.001), such that the cells can¬ 
not rearrange adequately. B The membrane force 
constant is relatively high ( F LB = 0.01), such that 
the rearrangement activity alters the topology sig¬ 
nificantly. C The membrane force constant is even 
higher ( F LB = 0.02). D The tissue shown in B is 
relaxed towards equilibrium. The color code denotes 
the domain identifier value (red = young cells, blue 
= old cells) 




Figure 13: Comparison of the Cell Topologies 
to Literature Results. P n denotes the occur¬ 
rence fraction of an n-sided polygon. Polygons with 
n > 9 are lumped into n = 9. A ’Low/medium/high 
tension’ correspond to different membrane tensions 
(. F LB = {0.001,0.01,0.02} in LB units), ’high ten¬ 
sion’ has a slightly higher occurrence of penta- and 
hexagons. B ’Farhadifar experiment’ is the polygon 
distribution in a third instar Drosophila wing disc. 
’Farhadifar case I’ and ’Farhadifar case II’ denote 
simulations with low and high relative contractility, 
respectively. ’Medium tension relaxed’ corresponds 
to an equilibrated simulation with F LB =0.01 after 
proliferation has stopped. 


23 












